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Abstract 



We perform a numerical study of the critical regime for the general relativistic collapse 
of collisionless matter in spherical symmetry. The evolution of the matter is given by the 
Vlasov equation (or Boltzmann equation) and the geometry by Einstein's equations. This 
system of coupled differential equations is solved using a particle-mesh (PM) method. 
This method approximates the distribution function which describes the matter in phase 
space with a set of particles moving along the characteristics of the Vlasov equation. The 
individual particles are allowed to have angular momentum different from zero but the 
total angular momentum has to be zero to retain spherical symmetry. 

In accord wih previous work by Rein, Rendall and Schaeffer, our results give some 
indications that the critical behaivour in this model is of Type I (the smallest black hole 
in each family has a finite mass). For the families of initial data that we have studied it 
seems that in the critical regime the solution is a static spacetime with non-zero radial 
momentum for the individual particles. We have also found evidence for scaling laws for 
the time that the critical solutions spend in the critical regime. 
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Chapter 1 



Introduction 
1.1 Gravitation and Critical Phenomena 

Einstein's theory of gravitation connects the geometry of spacetime to its matter content 
via the field equations^: 

G ab = 8nT ab . (1.1) 

Here G ab is the Einstein tensor, whose coefficients are complicated functions of the metric, 
g a b, and its first and second derivatives. T ab is the stress energy tensor, which depends 
on the matter, and also, in general, on the metric. 

In a given coordinate system, these equations are non-linear partial differential equa- 
tions for the metric coefficients and can give rise to very interesting solutions for the 
spacetime. Some of the most interesting phenomena occur at the threshold of black hole 
formation as was first discovered numerically by Choptuik who studied the collapse 
of a massless scalar field in spherical symmetry. 

For initial data characterized by one parameter p (the maximum amplitude of the 
scalar field for example), Choptuik found a critical value p* such that for values of p, 
p > p*, the evolution gave rise to the formation of a black hole, while for p < p* the scalar 
field dispersed to infinity. Near p = p* the space-time approached a universal solution, 
independent of the particular choice of the initial shape of the pulse (e.g. gaussian, 

tanh, etc). It was found that the near-critical dynamics was characterized by self-similar 

1 We will use units where c = 1 and G = 1 (c being the speed of light, and G Newton's gravitational 
constant) 
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oscillations with scaling factor e A . 

Choptuik also found that the masses of the black hole formed obeyed a scaling law: 

M BH <x\p-p*\i (1.2) 

where 7 was a universal exponent, again, universal in the sense of being independent of 
the choice of the family of initial data. This equation shows that, in principle, one can 
form black holes with arbitrarily small masses in the model. The transition to black hole 
formation in this case is thus continuous in the mass of the black hole. In analogy with 
the critical behavior in statistical mechanics this was called a Type II critical solution. 

A lot of work has been done in the last few years trying to find the critical solutions 
for different types of matter. Type I transitions (i.e. transitions where the mass of the 
smaller black hole is finite) have been also found |2[ . In this case, the solutions that have 
been found are either static or periodic ||. Here it was found that the time that the 
solution stays in the critical regime scales like: 

r ~ — a In \p — p*\ (1.3) 

where, again, a is a universal constant^} 

The process of tuning p to p* can be understood as choosing initial data such that at 
the threshold of black hole formation, we tune-out any growing component. Since we can 
tune away the growing components using only one parameter, the solution apparently 
has exactly one unstable mode. Linear perturbations of the critical solutions [§] agree 
with that observation. 

In this thesis we study the critical regime for collisionless matter in spherical symme- 
try. If we think of the matter as being composed of individual particles, allowing these 
particles to have angular momentum different from zero (keeping the net angular momen- 
tum equal to zero to remain in spherical symmetry) can give us some insight concerning 
2 Note that this constant depends on the units used. 
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the importance of angular momentum in critical phenomena, without needing to go to 
spacetimes with less symmetry. 

1.2 Collisionless Matter (A little bit of history) 

Collisionless matter (or dust) has been widely studied in general relativity. Einstein 
0] used this type of matter to investigate the physical significance of the singularity 
at r = 2 M for the Schwarzschild solution. He studied clusters of particles rotating in 
circular orbits (and spherical symmetry) and proved that for clusters in equilibrium the 
Schwarzschild mass function cannot reach M(r) = r/2 (r being the areal coordinate), 
therefore these configurations cannot have an event horizon nor a singularity. He postu- 
lated that this result also would hold for systems not in equilibrium and concluded that 
in "physical reality" there is no configuration of matter such that M(r) = r/2 (i.e. black 
holes cannot form). Although he came to the wrong conclusion, he studied one of the 
first models for relativistic stars. 

In the case when the angular momentum of each particle is zero, and the spacetime is 
spherically symmetric, Oppenheimer and Snyder || showed that a homogeneous distri- 
bution of this kind of matter generally forms a singularity. Using comoving coordinates 
they found the solution for the spacetime and deduced that singularity formation would 
happen in finite proper time (at least from the point of view of the comoving observers). 
Tolman || and Bondi generalized this result for non-homogeneous configurations, 
restricting to the case where no particle overtakes any other during the evolution. 

Eardley and Smarr || considered the same model as Tolman and Bondi and studied 
the existence of maximal (K = 0) and constant mean curvature (K = K D ) slicing condi- 
tions. Shapiro and Teukolsky ||, |T(|, |H| , |12| , |13| have studied this system in detail 



in the context of models of stellar clusters. They developed a code very similar to ours 
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to study properties of relativistic clusters, such as stability and collapse to a black hole. 
They also extended the code to treat more general spacetimes , i.e. with axial symmetry, 
|14| to study naked singularities. 

Rendall ]TjJ, |17| , and Rein fl5j have a series of papers addressing general properties 
for the Einstein- Vlasov system (the Cauchy problem, existence of static solutions in 
spherical symmetry, existence of general asymptotically flat solutions, etc.). 

More recently, Rendall, Rein and Schaeffer JTSJ have published the first study of 
critical behavior (black-hole threshold behavior) for this type of matter. Indeed, in this 
last paper Rendall et al argue that the critical solution in the case of collisionless matter 
is generically Type I. Our results agree with theirs in that we also find evidence for 
Type I transitions for certain initial data families which we have studied in this thesis. 
In addition to that result, we argue here that the critical solutions in this model are 



static, and we present evidence for scaling laws of the form (|1.3| ). Finally we have 
preliminary investigations aimed at determining whether the critical solution is unique 
(universality). 



Chapter 2 
The Equations 



The dynamical state of dust can be described by a distribution function, /: 

f(x a , Pa ) =dN/dV p (2.1) 

where N is the number of particles and V p is the volume in phase space. In our case, 
and since the matter is collisionless, the volume in phase space is a conserved quantity 
during the evolution of the system (Liouville theorem). This implies that the distribution 
function is also a conserved quantity: 

f - = (2.2) 
dr 

This is the collisionless Boltzmann or Vlasov equation. This equation plus Einstein's 
equations (|1 . 1|) , all restricted to spherical symmetry, i.e.: 

f(t,x i ,p i )= f{t,Rx\Rpi) with ReSO(3) i = 1,2,3 (2.3) 

form the system of equations that we want to solve numerically. 

Numerical solutions of PDEs generally involves discretization of the continuum prob- 
lem, and here we can consider at least two approaches: 1) Consider a set of particles 
which approximates the distribution function at the initial time and evolve the particles 



along the characteristics of equation (|2.2| ). The geometry is computed at each time step 
using the energy densities derived from the positions and velocities of the particles. 2) 
Discretize the Vlasov equation for the distribution function in phase space and solve for 
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it directly in phase space. In this case the energy densities, required in the update of the 
geometry, are calculated by integrating the distribution function over momentum space. 

In the particle approach, we have the problem that a given set of particles is just one of 
the infinite number of possible realizations for given density and velocity profiles, therefore 
we introduce statistical errors. In the second case we have to solve the Vlasov equation, a 
partial differential equation in phase space. Although we have taken the particle approach 
due to the simplicity of the evolution equations, we include the equations for the direct 
integration in Appendix B. 

As with any problem in "3+1" (or ADM) numerical relativity, we want to be able to 
specify initial data on a spatial hypersurface and then evolve these data in time. To do 
this, we need to split the equations Ql.lj) into a set of constraint equations (equations that 
must be satisfied at each instant of time) and dynamical or evolution equations (equations 
that tell us how to evolve the geometric quantities in time). In the most general case, we 
will have four constraint and six second-order-in-time evolution equations. 

In order to do this splitting, we will make use of the 3 + 1 formalism due to Arnowitt, 
Deser and Misner (ADM), for a review of this formalism see PD| . We have followed a 
summary by Choptuik [5I| which is itself based on a review by York p4| . 

In the remainder of this chapter, we derive the system of equations that we solve. We 
start with a summary of the 3 + 1 formalism and its restriction to the spherical symmetric 
case, and a definition and description of the maximal-areal coordinate system. In section 
2.2 we explain how the Einstein's equations are coupled to the matter, and in 2.3 we 
explain the particle evolution equations in our coordinates. The chapter ends with some 
calculations of conserved quantities. 
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2.1 The 3+1 Formalism of General Relativity. 

We consider a spacetime manifold with metric g a f, defined on it, and assume that we can 
slice the manifold into spacelike hypersurfaces, defined by t = constant, where t is our 
time coordinate. At least locally then, the space-like surfaces can be described in terms 
of a closed one-form, Qj^, which is just the gradient of t 

n a = dt = V a t. (2.4) 

The norm of this one-form is: 

g ab n a n b = -a~ 2 } (2.5) 

where a is a scalar function commonly known as the lapse function. As long as our slices 
are spacelike, a will be strictly positive. We introduce a normalized one-form, n a , defined 
by 

n a = -aQ a = -aV a t. (2.6) 

The contravariant vector, n a = g ab n , which is the future-directed normal to the surfaces, 
can be viewed as the four velocity field of observers moving orthogonally to the slices. 

In this section, again following York, we will describe the 3+1 split in terms of 4- 
dimensional tensors. At the same time, we will need to decompose various tensors (in- 
cluding the Einstein and stress energy tensors appearing in fll.lQ ) into pieces parallel to 
n a ("timelike" pieces) and pieces orthogonal to n a ("spacelike" pieces). Thus we define 

W n = _ W a n ^ (2.7) 
where the sign convention is again York's, while for covariant vectors we define 



W n = W a n a . 
Indexes in this section are abstract indexes as in [E6| 



(2.8) 
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For the projection onto the hypersurfaces we use the projection operator: 

r b = S\ + n a n b (2.9) 
The projection of the space-time metric to the hypersurfaces is: 

lab = l C al d bgcd = 9ab + n a 7l b (2.10) 

and is, in fact, the metric induced on the hypersurfaces by the 3+1 splitting. (Note that 
7 ab is not necessarily the inverse of j ab , as follows immediately from the definition of 7%, 
and in this section we raise and lower all indexes with the four dimensional metric g ab 
and its inverse). We can now introduce the natural derivative operator on the space-like 
surfaces, by projection of the four dimensional covariant derivative: 

D a = 1 b a W b . (2.11) 

It is easily shown that this derivative is compatible with the spatial metric: 

D albc = 0. (2.12) 

Having defined a derivative operator on the hypersurfaces, we can compute the asso- 
ciated Riemann tensor 3 7Z abc d which measures the intrinsic curvature of the hypersurface. 
For example, for a spatial one form u a {u) a n a = 0) we have 

(D a D b - D b D a ) u c = 3 n abc d iu d . (2.13) 

We can then also define the spatial Ricci tensor, 3 1Z ac = 3 TZ abc b and the spatial Ricci 
scalar, 3 TZ = z Tt a . 

Finally, we define the extrinsic curvature tensor, K ab : 

K ab = -D a n b = V a n b - n a a b , (2.14) 
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where a b = n a D a n b is the dual to the four acceleration field of the observers moving with 
the slicing. 

The geometry of the space-time is completely defined by the spatial metric, r y a b, which 
describes the geometry of each slice, and the extrinsic curvature tensor which tells us 
how each slice is embedded in the four dimensional space-time. 

We can write Einstein's equations in terms of the tensors defined above, as well as 
the following projections of the stress energy tensor: 

p = T ab n a n b (2.15) 

f = 7 vtV) ( 2 - 16 ) 

gab = Y^b^cd ( 2 _ 17 ) 

These are interpreted as the local energy density, momentum density and spatial stress 
tensor, respectively, for an observer moving orthogonally to the space-like hypersurfaces. 
The Hamiltonian constraint is found by contracting Einstein's equations with n a twice, 
G ab n a n b = T ab n a n\ yielding 

3 TZ + K 2 - K a b K\ = lQnp (2.18) 

where K = g ab K ab = r y ab K ab is the trace of the extrinsic curvature and 3 1Z is the spatial 
Ricci scalar defined previously. 

If we contract Einstein's equations once with n a , and then project the resulting vector 
onto the hypersurface, 7 a c G cb nfe = 7 a c T cfe nb, we get the momentum constraint: 

D b K ab - D a K = 8irj a . (2.19) 

These two constraint equations involve only spatial tensors and spatial derivatives of 
spatial tensors, and must be satisfied on each spacelike slice, including the initial slice. 
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In order to derive the evolution equations we use the Lie derivative along the vector 
field, N a : 

N a = an a + f3 a = d/dt (2.20) 

where j3 a is an arbitrary spatial vector field, commonly known as the shift vector. Note 
that N a satisfies: 

N a tt a = 1 (2.21) 

It can be shown that Lie differentiation along N a commutes with the projection operator. 
We can now write the rest of equations (|1.1|) as two different sets: (1) the evolution of 
the spatial metric: 

C aab = -2aK ab + C plab (2.22) 

which can also be viewed as a definition of K ab , and (2) the evolution equations for the 
extrinsic curvature: 

L t K\ = C p K a b - D a D b a + a 3 K a b + KK\ + 8tt {^P\ (S - p) - S a b ^j (2.23) 
where S is the trace of tensor given by equation Q2.17|) . 

2.1.1 3+1 Formalism in Spherical Symmetry 

We now restrict our attention to spherical symmetry. In contrast to the previous section, 
where our tensor expressions involved abstract 4-dimensional indices, we will be mostly 
concerned with the components of specific 3-tensors in this section. Thus, Latin indices 

k, • ■ ■ range over the spatial values 1, 2 and 3, and all indices of the spatial tensors 
are raised and lowered with 7^ and 7^ respectively, where 7* J 7-,-fc = 5\. 

The most general 3-metric in spherical symmetry can be written as: 

V = a 2 (t, r) dr 2 + r 2 b 2 (t, r) d£l 2 (2.24) 
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where dfl 2 is the metric on the unit sphere, dVt 2 = d9 2 + sin 2 9d(p 2 . From the 3-metric, 
we can compute the associated connection coefficients or Christoffel symbols: 



r k ~ 2^ (^ n 3< k 7nfc,jr' ljk,n) 



(2.25) 



In our coordinates, the non-zero connection coefficients are: 



V\ r = a! I a Veo = ~ (W) / (2a 2 ) T r H = - sin 2 9 (r 2 b 2 )' / (2a 2 ) 

r e r9 = ((r 2 6 2 )') / (2r 2 6 2 ) V s = -sm9 cos 9 



r% = (r 2 6 2 )' / (2r 2 6 2 ) 



T^m = cot 9 



We also need to compute the components of the Ricci tensor: 



(2.26) 



3p _ pn p« I pn pm pn pm 

*Hj L ij,n L in,j ~r J- nm J- ij J- jm J- in 



(2.27) 



from which we find the non-zero components: 



(r 2 6 2 )' / (r 2 6 2 )l ' + a' {r 2 b 2 )' / {ar 2 b 2 ) - \ \{r 2 b 2 )' / (r 2 6 2 ) 



R ee = - (r 2 b 2 )' I (2a 2 ) • + 1 - a! (r 2 6 2 )' / (2a 3 ) 



(2.28) 



i R thth = sin 2 9R 



The mixed components are given by 



3 R\ = -2 



{rb)' I a I (arb) 



R\ = R<t> 4> = a - (rb {rb)' /a)' / (a {rb) 2 ) . 



(2.29) 
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Another quantity which we need is a^, where "|" denotes covariant differentiation with 
respect to the 3-metric. 

a|i |j ' = y fc (a )ifc -r n ifc a Jn ) (2.30) 
which has the following non-zero components: 

a{r \r = ( a ' /a)' /a a,/ = a' (r 2 b 2 )' / (2r 2 a 2 b 2 ) a,/ = (rb)' a'/ (rba 2 ) (2.31) 
Finally, the Ricci scalar, 3 R = ^R^, is given by: 



arb 



_ 2 ({r^_\ Jr a L l 1 (rbf 



(2.32) 



a ) ' rb \ a 2 

At this point we have calculated all the geometric objects intrinsic to the spatial slices 
that we need. Now, we embed these slices into the most general spherically symmetric 
spacetime whose metric can be written: 

ds 2 = (-a 2 + a 2 !3 2 ) dt 2 + 2a 2 (3dtdr + a 2 dr 2 + r 2 b 2 dVi 2 (2.33) 

where a, (3, a and b are all functions of r and t. Note that due to the spherical symmetry, 
the shift vector has only a radial component; f3 l = (/?, 0, 0), so (3 will be called the shift 
function. In these coordinates n M = (—a, 0, 0, 0) and n M = (1/ct, — (3 /a, 0, 0). 
The extrinsic curvature can be written in terms of the 4-metric coefficients: 

Ka = ~ - Hj,k(3 k - l ik (3\ 3 - lkj (3 k )j (2.34) 

From this we compute the non-vanishing components for the current spherically sym- 
metric case: 

K„. = ?-({a(3)'-a) (2.35) 
K e6 = r -^(p{rb)' -rb) (2.36) 
K H = sin 2 6K 0e (2.37) 
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The corresponding components with mixed indexes are: 

\ap)' - a) (2.38) 













aa 




1 


K% 






rba 




= K% 



(rb) 1 - rb) (2.39) 

(2.40) 



and the trace of the extrinsic curvature is 



*- = (MM + 1(^-1) ™ 

Finally, the only non-zero component of DiK = K\i is 

K V = K r = - — (- + -) {ap - a) + — {{a(3f - a') - \ol U^- - l) 
aa \ a a J aa a 2 \ rb J 

2.1.2 Maximal Areal Coordinate System. 

Einstein's equations allow for coordinate freedom that we need to fix. We have chosen 
maximal-areal coordinates (defined below), but have also included the equations for the 
polar-areal system in Appendix C. The main advantage of maximal-areal coordinates is 
that, in contrast to the polar-areal case, the slices used can penetrate apparent horizons. 
As the name suggests, in the maximal-areal coordinate system, the radial coordinate is 
areal, so that the proper area of 2-spheres with radius r is Attt 2 . In terms of the general 
spherically-symmetric 4- metric (|2.33| ), this means that b(r,t) = 1. The time coordinate 
is fixed by demanding that the 3-slices be maximal, i.e. that K (r,t) = 0. This leads to 
a condition on the lapse function a (r,t) (the so-called slicing condition), which must be 
satisfied at each instant of time. 



Thus, in maximal-areal coordinates the 4-metric ( [2.331 ) takes the form: 



ds 2 = (-a 2 + a 2 (3 2 ) dt 2 + 2a 2 pdtdr + a 2 dr 2 + r 2 dtt 2 (2.42) 
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Using the expressions we derived in the previous section, we compute the relevant geo- 
metrical quantities we need in order to write Einstein's equations in terms of these metric 
coefficients: 



R = 2 



-2 (1/a)' + a/r (a - l/a 2 )l /(or) 



R\ = -2(l/a)'/(ar) 



= ify = [ a - (r/a)'l /(or 2 ) 



K r r = ((a/3)'-d) /(a a) 



K e e = K^ = /3/ (ro) 



(2.43) 



ifVi; = K r rjr + 2 - /r 



l£| r = 



a| r ' r = (oi'/a) I a 

We can now specialize equations ( |1.1| ): 
Hamiltonian Constraint: 



a/ = a^ = a'/(ra 2 ) 



— = \ a 2 rK e e 2 + 4vrra 2 p + (l 
a 2 yr \ 



Momentum Constraint: 



Evolution of the 3-metric: 



2r 



r 



(2.44) 



(2.45) 



a = 2aaK% + (a/3)' 
P = arK\ 



(2.46) 
(2.47) 



Slicing Condition: 



a" = a' (- --) + ^ (2r- + a 2 - l) + 4vra 2 a {S - 3p) 



a r / r° \ a 



(2.48) 
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where in deriving the last formula, we have made use of the slicing condition (K = 0), 
and the evolution equations for the extrinsic curvature components, equations (|2.23| ). As 
we will discuss in more detail below, we have chosen to do a constrained evolution, which, 
in this case means that we use the constraint equations, rather than evolution equations, 
to update a and K 6 e . 

2.2 Stress-Energy Tensor 

In this section we explain how we calculate the energy densities p, j r and S that appear 
in equations ( [2.44| - |2~4"8|) . We approximate the distribution function (|2.1|) by a set of N 
"spherical particles". Since the particles only interact with each other gravitationally, we 
have 

N 

T» v =Y, T t ( 2 -49) 

i=i 

where Tf v is the stress energy tensor for a single particle. For a point particle we can 
write: 

Tr = ^-S(r-fi(t)). (2.50) 

mi 

Here, pf is the p component of the 4 momentum of the z-th particle, m, is its rest mass, 
and fi(t) is the spatial position of the particle at time t. 

Using equations ( |2.15| - |2.l7| ) for p, S and j r specialized to our coordinates, we get: 



[p\ = a 2 [T tt \ (2.51) 

b"r] 4 = (2-53) 

Now we relax the point particle approximation and assume that each particle is a spher- 
ically symmetric shell of mass, uniformly distributed over a region Ar of spaceQ. These 

2 Later on this Ar will be the mesh spacing used in the finite difference solution of the geometrical 
equations. 
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shells of matter are averages of the shells at that point with magnitude of the angular 
momentum I pointing in all posible directions. Therefore the angular momentum of the 
average is zero, I = 0, but I 2 ^ 0. The proper volume that each particle occupies is then: 

t 2 t 

Vi = — Ar \ J^g dcfidO = AnAraa^-^ (2.54) 

mi J rrii 

and we can approximate the delta function that appears in ( 2.50Q by l/V*. This yields: 

»< = sLf (2 ' 55) 

St = [S%] ( + [S"J, = 4;rAro 3 [p,J ( ' r 2 + ^Ararflp'l (2 ' 56) 

Here a, a and /3 are evaluated at r = r«, [p*^ is defined as [p t ] i = a \p\, and [l]^ = 
[pe]j 2 + [p^]// sin 2 is the square of the magnitude of the angular momentum of the z-th 
particle. We then introduce quantities which do not explicitly depend on the geometry 
(since after updating the particle positions we want to solve for the geometry): [p],. = 
a \p] it [S r r ]. = a 3 [S%] 4 , [S\] . = a [£»J. and fr]. = a \j r ] t . 

We interpolate the one-particle quantities to the continuum and sum over all the 
particles to find the total values: 

N p 

f = ifiW(r-n) (2.58) 
i=i 

where fi is any of the barred quantities defined above for each particle, / is the corre- 
sponding quantity in the continuum case, and W(r — rj) is an interpolation function (see 
section fj.4| for further explanation of how we define this function). Having defined Q2.58| ) 
we can now write equations ( |2.44| - |2.48"D as 

1 - " 2 + lra 2 K e e 2 + Aixarp (2.59) 



a 2r 2 
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K%' = -Ir\-^ (2.60) 
r a 



6 ' 




e — 


A g 




r 








,(a> 


a" = 


a 




y a 


P = 


arK e 6 



+ ^ ( a 2 - 1 + 2r- ] + 4vraa I ^ + S a a - 3p ) (2.61) 

\ (XI \ (X 



(2.62) 



2.3 Evolution Equations 



The equations of motion for the particles are just the geodesic equations (the character- 
istics of the Vlasov equation). These can be expressed in terms of the four momentum 
of the particle: 

p a \7 aP b = 0, (2.63) 

so in our coordinate system we have: 

^ = -T\ t (p t ) 2 -2T t tr fp t -T t rr (f) 2 -T t ee ^ (2.64) 

-nr I t\ 2 orr rt rr / r\2 rr ^ 



^- = -r tt (/) -2rvpV-r r rr (pT-r r w^ (2.65) 
^ = -2r%pV-rW (2.66) 

GST 

where the T's are the Christoffel symbols defined by ( |2.25|) , and I 2 = Pe 2 +P<j> 2 / sin 2 6 is the 
square of the magnitude of the angular momentum as before. We want to consider the 4- 
momentum as a function of our time coordinate and not the proper time of each particle. 
In order to do the transformation we can use the chain rule [djdr = (dt/dr) (d/dt) and 
p l = dt/dr): 



% = -V\ t p t -2Y\ T f-V\ r ^-V t ee 4j (2.67) 
dt p l p l r^ 



dn r (n r ) 2 J 2 

-j- = -V r aP t -2V r tr p r -V r rr ^--V r ee—^ (2.68) 
dt p l p'r 4 

IF ~ - 2T ^—r- T ^-r ( 2 - 69 ) 
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Furthermore it is convenient to express these equations in terms of p r rather than p r : 

Pr = grtff = a 2 (t, r)/3(t, r)p l + a\t, r)p r (2.70) 

Solving for p r : 

p r w=sB r )~ i3{t,r)pt{t) (2 ' n) 

We can now compute the total derivatives with respect to coordinate time as follows: 

d d dr dr d d p r d 

— = 1 = h (2.72) 

dt dt dr dt dr dt p % dr 



Applying this operator to equation (|2.71 ) we get: 



df_ = ]_dPr_ 2 f_ + _pdp__ p t (^E + r ^_^E\ (273) 

dt a? dt a 3 \dt p t dr J dt \dt p t dr J 

Substituting equations ( |2.71| ) and ( |2.73| ) into equation ( |2.68| ) we obtain: 

d JL = «V + f Vr + 1*^ + JL (2.74) 
at dr dr cr dr p t p l r A 

which is the evolution equation for p r . To derive the evolution equation for r we use the 
definition of p r (p r = dr/dr) which after some manipulation yields: 

- = ^--/3 (2.75) 
dt a 2 p t 

To find the time component of the 4 momentum we make use of the normalization 
condition p^p^ = —m 2 : 



ap t = Jm 2 + ^ + ^ (2.76) 



It is also convenient, as previously mentioned, to use p* = ap l rather than p l itself. Then 
the resulting geodesic equations in these coordinates are (we have included equations for 
p^, p^ and dp^ /dt because use we use them sometimes for visualization purposes): 

dp r da_ t d/3 adap r 2 I 2 a 

—rr = -JTP + ~a~Pr + -^7T~T + ( 2 - 77 ) 

dt dr dr cr dr p l p t r i 
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p> = ( / m2 + ^ + il (2 . 79) 



d0 op* 



(2.81) 



/,2 



1 // 2 



= —n \-~P ( 2 - 82 ) 

sm v \ / 



2.4 Conserved Quantities 



(2.83) 



Here we include the calculation of two quantities that we have used to check the code. 
The first quantity that we discuss is a conserved quantity in flat spacetime related to the 
energy of the particles. This gives us a check of the evolution when we fix the geometry 
to be Minkowskian (flat), which was useful in the development of the code. The second 
quantity that we calculate is a mass aspect function that coincides with the ADM mass 
at infinity. The conservation of the value of this function at the outer limit of our range 
of integration gives another check in the case when we solve the fully coupled problem. 

In flat space-time, the geometry is static, and therefore we have a time-like Killing 
vector field. Choosing coordinates such that: 

ds 2 = -dt 2 + dr 2 + r 2 d6 2 + r 2 sm 2 6d<f) 2 , (2.84) 

we can write this Killing vector field as = d/dt. If we contract with the stress 
energy tensor T^ u , we get a vector field whose divergence is zero: 

V M (T^%) = (V^TH + (V^,) = 0. (2.85) 
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The first term on the right hand side is zero by the conservation of the stress-energy 
tensor, and the second is zero by the assumption that satisfies the Killing equation^, 
and by the symmetry of T^ u . Integrating this divergence over the four volume we get: 

/ V,(Ta)^j;= / T^^^~gd 3 x = (2.86) 

JV JdV 

where the second integral is an integral over the 3-hypersurface dV (boundary of V). 
We choose this 3-hypersurface to be composed of two different hypersurfaces of constant 
t, t = t\ and t = t2, joined by a timelike-hypersurface at spatial infinity. We assume 
that the stress-energy tensor vanishes at spatial infinity, so that it does not contribute 
to the integral. Let us also assume that = (—1,0,0,0) for the surface with ti, and 
n v = (1, 0, 0, 0) for t 2 . In this case: 

/ T tt J(^d 3 x= ( T tt \[Wgd?x (2.87) 

JdV H JdVt 2 V 

and since t\ and t 2 are arbitrary this means that the integral is constant for any value of 
the time and the integrand calculated using ( |2.50| ) 

N v t2 

(2.88) 

is a conserved quantity 0. 

The other quantity that we calculate in this section is the mass aspect function: 

M = 5( 1 + 54)=5( 1 + W -?l (2 - 89) 



To derive ( [2.89 ), we consider the transformation relating the metric (|2.42 ) to the Schwarzschild 



metric (in a region of vacuum): 



1 - ^) dt 2 s + (l - ^) 1 dr\ + r s W (2.90) 



3 V M £„ + V„£ M = 



4 Not only ( p. 88 ) is a conserved quantity but, in flat spacetime, p* for each individual particle is 
conserved. 
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Here the V subscript stands for Schwarzschild and we wish to stress the point that the 
Schwarzschild time and our time coordinate are different, while the radial coordinates 
are the same. Following || we get: 

dt s = C(dt + Ddr) (2.91) 

dr s = dr (2.92) 

where the last equation holds because both radial coordinates are areal. Using these 
formulas in ( |2.90p we get: 

( 1 - —) C 2 dt 2 -2(l-—) C 2 Ddtdr 
r J V r / 

1 - —J C 2 D 2 dr 2 + ( 1 - —J dr 2 + r 2 dfl 2 (2.93) 
It is easy to see that in order for expressions ( p.42j ) and (|2.93|) to be equal we need: 

C 2 = ~ ~fXff (2-94) 
1 - 2M/r v ; 

D = fl, R 2 ( 2 - 95 ) 
— cr + a z p z 

If we now compare the dr 2 terms of ( |2.42| ) and ( [2.93D , we get an equation for the mass: 

" 11 + 1 = « 2 > ( 2 - 96 ) 



-a 2 + a 2 (3 2 \-2M/r 

which after some manipulation can be rewritten as: 

2M 1 (3 2 , ms 

— = 1 - - + - 2 • 2.97 

This last equation is interesting in its own, because when this quantity, 2M(t,r)/r, 
becomes equal to 1, we know that a marginally trapped surface has been formed. We 
can see this by computing the expansion of the null geodesies (see for instance p5|), 
which in this coordinates can be written as: 

1 - a(t, r) r K° e {t, r) = 1 - a ^f^ r \ (2.98) 

a{t, r) 

Therefore, if the expansion is zero, 1/a 2 = (3 2 /a 2 , and 2M(t,r)/r = 1. 
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The Code 



In this chapter we explain the details of the method we have used to solve the equations 
which were presented in the last chapter. 

3.1 Particle-Mesh Methods 

We have used a method of calculation known as "particle-mesh" (PM). For a very com- 
plete review of this and other particle methods see [^T[ . 



This technique solves the coupled Einstein- Vlasov equations by splitting each time 
step into two stages: 1) solution of the field equations on a mesh, and 2) updating of the 
positions of the particles via their equations of motion. 

To provide some motivation for the use of PM methods let us first consider the iV body 
problem in Newtonian gravity (where we assume iV is very large). We could calculate 
the forces experienced by each particle Fj directly using 

Once the acceleration at each particle is known, we could integrate the equations of 
motion over some small time interval At giving the new positions and velocities of the 
particles. This method is denoted particle-particle, and computationally is very expensive 
since we have to perform 0(N 2 ) calculations per time step. 

The "particle-mesh" (PM) approach to the problem involves calculation of the mass 
density due to the N particles on a mesh (assuming each particle occupies a finite region 

22 
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in space). On this mesh, we can use finite difference approximation (or some other 
discretization technique, like a spectral method) to solve the Poisson equation for the 
gravitational potential: 



Once the potential is found, we can interpolate the value of its derivative to the actual 
particle positions to find the forces. A time step is then concluded by solving the evolution 
equations for each particle to provide updated positions and velocities. 

PM methods are most effective when the close interactions between the particles are 
not important for the evolution of the system. The main reason for this is that smoothing 
of the particles in a finite region of the space produce a large error in the density close to 
the particle positions. However, in cases when the dynamics is dominated by a mean field 
the particle-mesh method is much faster than direct integrations like particle-particle. 

In our case we want to solve the Vlasov equation (equation (|2.2j) ) by approximating 
the distribution function with a set of point particles. Further approximation is necessary 
to couple to the geometric equations, since a point particle gives rise to infinite densities. 
To solve the geometric equations we thus have to find the stress-energy tensor for the 
particles, where we consider each particle to be smoothed out over a finite region of space. 
That is precisely the PM method defined above. 

In the next sections of this chapter we will explain the specifics of the integration of 
the field equations, the evolution equations for each particle and the interpolation scheme 
that we have used. We also include pseudo-code for the main routine of our program in 
Appendix D. 
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3.2 The Field Equations 



We first explain how the equations ( |2.59[]2l)2j ) for the geometry are solved numerically, 



assuming we know the quantities p, j r , S^., Sg []. The first two equations, equations 
( P-59| - |2~]60| ), are integrated from the origin, r = 0, using the LSODA |22j integrator. The 
boundary conditions are given by the spherical symmetry of the space-time, and by 
the demand that the spacetime be locally flat at r = 0. They are a(t,0) = 1, and 
K e e (t,0) = 0. 

We compute the values of the functions and K e oi on a uniform grid of N r points 
at positions rf r = 0, r 1 = h, r 2 = 2h,..., r Nr = r max . 

In order to compute the values at r = r i+1 , we supply to LSODA the values of the 
functions at r = r i and the derivatives computed using equations ( |2.59| - |2.60| ) at r = r i+1 / 2 , 
using p and j r averaged between the i-th point and the (i+l)-th point: 



[Pirn/2 = \ (\P]i + lP] i+ i) (3-3) 

b"r] i+1/2 = \ + (3.4) 



2 
1 

2 

Once we have calculated a we can solve the slicing equation: 



n' 2 In n' 
a" = a '( ) + -i(a 2 - 1 + 2r-) + 4naa(S r r + S a a - 3p), (3.5) 



a r r 2 a 



with the boundary conditions: 

a'{t,0) = (3.6) 
a(t, oo) = 1. (3.7) 



Here the first condition follows from ( j3.5|) demanding that a be regular at r = and 

the second one follows from asymptotic flatness, and the demand that t measure proper 
-'^we will explain how to compute these quantities in section 3.4 
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time at infinity. We use a second-order finite difference approximation on the previously 
defined mesh of points {r^}: 

i-1 Qfi+l — Oli-i ( Oi+x — di-x 



Qfj+i — 2aj + a 

h 2 " 



2A 



2a* 



1 + 2r. 



2/m; r. 



= 1 + 



+ 



2/ia< 



- 3 [4 



+ 



Rearranging this equation gives us: 



T? + 9i ) «i + f A - 77T J = 



/i 2 



/i 2 2/i 



where: 



ft 



Oj+l — 0>i-l 

2/ia,- 



2 

r,- 



a; — 1 + 2r, 



2/ia, 



47ra,- 



+ 



3 \P\ 



In addition to these equations we have the boundary equation at r = 0: 

a 2 = 



l//i 2 + / 2 /(2/i) \ 

+ v/i 2 -/ 2 /(2/i); 



/ -2//i 2 -^ 2 \ 
ai+ ( 4+ l/^-/ 2 /(2fc)J 



(3.8) 
(3.9) 

(3.10) 

(3.11) 

(3.12) 
(3.13) 

(3.14) 



which can be derived from the forward finite difference approximation to a' = at r = 0: 

— 3«i + 4a 2 — ct3 



2h 







and equation ( |3.11| ) with i = 2. We have also the boundary condition at 



r = r r , 



(3.15) 



(3.16) 



N r 



This is an approximation to (|3.7|) at a finite value of r = r max , where M is the mass aspect 
function defined by equation (|2.89|) . These equations form a linear system of algebraic 



equations that can be solved using a tridiagonal solver (we have used the LAPACK p3 
routine dgtsv). 
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3.3 The Evolution Equations 



To evolve the particles' positions and momenta we integrate the geodesic equations ( 2.77 



2.78|) . The values of the coefficients in these equations (basically products and quotients 
of a, a, j3, a', a' and (3') must be calculated at the particles' positions using the values 
obtained at the mesh points {r^}. The mesh values are interpolated to the particle 
positions using the same operator kerner used to produce mesh values from particle 
quantities (this procedure is explained in section |3^) . These equations are integrated 
using the LSODA integrator as in the case of the Hamiltonian and momentum constraints 
( |2.59| - |2"1)H| ). For a given particle at position r n and momentum p n r at time step n, we 



calculate the new position r n+1 and momentum at t = t n+1 by suplying to LSODA the 
values of r n , p n and their derivatives at time step t = t n . This procedure has an accuracy 
O(At) (since the derivatives are calculated using the metric coeficients at t = t n ) where 
At = t n+1 — t n . In our program we chose a value of At proportional to h, i.e. At = Xh, 
where usually A = 1.0. 

We need to take special care if a particle leaves the range of integration (r > r max ) or 
if it reaches the origin. In the first case we stop considering the particle and we change 
the total number of particles to (N — 1). When a particle p reaches the origin, in other 
words when r p < 0, we reflect the particle, i.e. we set: 

[r], =► -M p (3-17) 

[Pr] p =► -\Pr) p (3.18) 

[*] P =► P], (3-19) 



As mentioned previously, equations ( |2.82|) , for p^, and (|2.83|) , for 0, are also sometimes 



solved for visualization purposes and is also integrated in the same way. 
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3.4 Interpolation and Restriction 



In this section we explain how we calculate the energy distributions on the mesh from 
a given set of particles. In other words we will discuss the different possibilities for the 



interpolation operator W(r—ri) appearing in equation ( 2.58 ). At the end of the section we 
explain how to use the same operators to restrict the coefficients of the geodesic equations 
calculated on the mesh (i.e. the geometric quantities), to the particles' positions. 

We can define a hierarchy of interpolation operators which can be classified by how 
many derivatives of the resulting density (interpolant) are continuous. The lowest-order 
interpolation (the resulting density is discontinuous, piecewise constant), called nearest 
grid point interpolation, assigns the density of each particle to its nearest grid point 
(NGP). Thus the interpolation kernel can be written as: 

1 : \x — x v \ < h/2 

P (3.20) 

: otherwise 

where x p is the position of the particle and h is the mesh spacing of the grid on which 



W(x - x p ) 



W(x-p) 







h/2 





i-1 



i+1 



Figure 3.1: Kernel for the first order interpolation (NGP). 

we solve the field equations. 

Figure |3.1| displays a typical NGP kernel operator for a particle at position p, and 
a grid with grid points at r = r^x, rj, r i+1 , etc. This kernel will produce a density 
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distribution given by: 



N 



p(?) = J2pp w ( x - 



X r 



(3.21) 



i=l 



A commonly used higher order interpolation is called "cloud in cell" (CIC) interpo- 
lation which distributes the density of each particle into two grid cells. In this case the 
resulting density is piecewise linear and the kernel takes the form: 

l — \x — x p \/h : \x — x p \ < h 
: otherwise 



W(x — x p ) 



(3.22) 



Figure pO is a graph of the kernel of the CIC interpolation. In our calculations we have 




i-l 



i+l 



Figure 3.2: Kernel for the second order interpolation (CIC). 

used the CIC interpolation operator. 

To restrict the Christoffel symbols (and other geometric quantities) calculated on the 
mesh to the particles' positions we use the same kernel in the following way: 



N r 



F(x P ) = F(xi)W(xi - x p ) 



(3.23) 



i=i 



where x p is the position of the particle, Xi are the grid points, and F is any of the 
coefficients that appear in the evolution equations. The coefficients F are products and 
quotients of the metric coefficients and their derivatives. In order to calculate derivatives 
we use the following finite difference approximation on the mesh: 



W\i = (ft+i - Qi-i) / h + (h 2 ) 



(3.24) 
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and then use equation (|3.23|) to find an approximate value of q' at the particle's position. 



3.5 Initial Set of particles. 

To initialize the set of particles which we evolve, we specify the particle distribution 
(number of particles per unit of areal coordinate) and the velocity distribution (specifi- 
cally the number of particles per p r and I) . This correspond to a distribution function in 
the sense of fl2.1|) which is separable in the r, p r and I variables. In other words: 

f(r,p r ,l)=R(r)P(p r )L(l) (3.25) 

Moreover, instead of specifying P as a function of p r we specify P = P(p r ) where 



Pr — p r / a - This allow us to calculate the value of p l = ym 2 + p 2 , + l 2 /r 2 , and therefore 
p, j r and S, for each particle without knowing the geometry in advance. This allows us 
to decouple the tasks of specify initial data for the particles, and ensuring that the initial 
data satisfies the constraint equations. 

We use 1-dimensional Monte Carlo techniques applied to each of the functions R(r), 
P(p r ), L(l) to get an specific ensemble of the particles. In theory the statistical error 
that we produce goes like 0(N^ 1 ^ 2 ), where N is the number of particles that we use to 
sample the distribution function (see section |4~3|). 



Chapter 4 



Checks 



In this chapter we explain different checks we have performed to test that we have properly 
implemented the algorithm. The first of these test was checking if the quantity given 
by equation (|2.88 ) was conserved in the flat spacetime case. In our code this value was 
conserved up to the tolerance of the LSODA integrator for initial data which was gaussian 
in both the radial and velocity distributions. We have also checked that a test particle 
moving in a fixed background follows the geodesic of the background, and that its energy 
is conserved. 

The rest of the checks that we describe here are for the fully coupled problem. In 
section |4.1| we discuss the conservation of the mass aspect function at the outer boundary 
(r = r max ). We also have simulated certain clusters of particles in circular motion known 



as Einstein clusters; results from these simulations are discussed in 4.2 in section 4.2 



At the end of the chapter we study how the numerical errors scale with the number of 
particles we use to sample the distribution function /, and as well as with the mesh 
spacing h. 



4.1 Mass Conservation 

The value of the mass function given by equation ( |2.89| ) at infinity should be a conserved 
quantity if the spacetime is asymptotically flat (see for instance |19|). As an approxi- 
mation to the value at infinity we can monitor the value at the maximum value of the 
radial coordinate, M(t, r max ). This check will make sense as long as no matter leaves the 

30 
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range of integration (note that energy cannot escape in the form of radiation since that 
is forbidden in spherical symmetry). In Figure |4.1| we plot M(t, r max ) for three different 
mesh scales (N r = 64, N r = 128 and N r = 256) and with a fixed number of particles, 
N = 10000. 



0.8 



0.6 

x 
E 

ii 

*j 0.4 



o a 







10 20 30 40 50 

t 

Figure 4.1: Value of 2M(t,r)/r at r = r max as a function of time. 

We observe that the value of M is conserved until t ^ 40 which s the time when the 
particles begin to leave the range of integration. In the inset we see a detail of the function 
between t — and t = 40 showing that the conservation is better than 1% for N r = 64 
and improves somewhat as we increase the number of grid points. Theoretically (and 
in the limit N — > oo) we should observe that the variation of the value of mass should 
decrease linearly with the number of points, i.e. first order. Although the convergence 
test here is not as definitive as it tends to be for pure finite-difference methods, the inset 
of Fig. fO| provides fairly good evidence that our code's mass conservation would converge 
linearly with h in the continuum limit. 
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■ N r =128 
' N r =256 
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4.2 Einstein Clusters 

Probably the strongest check of the reliability of our code is its ability to simulate clusters 
of particles that only have rotational motion (zero radial component of the 4- momentum). 
Let us calculate what initial conditions we must set in order to obtain such a cluster. 

Assuming we know what the energy density of the cluster of particles is, let us cal- 
culate the velocity profiles required in order to produce a cluster in equilibrium. In the 
process of obtaining the velocity profiles, we will also find the geometry of the space- 
time created by the cluster. The metric coefficient a(r) can be calculated using equation 
( 2.89| ). If the spacetime is static then (5{r) = and the resulting equation is: 



where we calculate M(r) by integrating the density p(r) over the radial coordinate, 
M(r) = / 47rr 2 p(r)dr. To calculate the lapse function we can make use of the evolution 
equation for the extrinsic curvature Q2.23 ), which in this case reduces to: 



— j = aa — 4-naar p. (4.2) 

After a little bit of manipulation this equation can be written as: 

l~a 2 _ M(r) 1 
a 2r r 2 (1 - 2M{r)/r)' 1 ' 1 

We want to find the angular momentum l(r) needed to get the particles rotating with zero 
radial component of the 4 momentum. Imposing that condition (p r = and dp r /dt = 0) 
on the geodesic equation for p r , equation ( |2.77| ), we obtain: 

da . l 2 a n .. A . 

~ TP + T=i = ( 4 - 4 ) 

Taking into account that 



P t = \ m 2 + l — , (4.5) 
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and using equation ( f4.3| ), we find the equation for l(r): 



2J1 



I (r) = m r 



'M(r) 



1 



(4.6) 



r J (1- 3M(r)/r) 

This is the same expression as the angular momentum for a particle in circular orbit 
around a Schwarzschild black hole with mass M(r) (see p6|). In Figure |4]^ we have 
plotted the areal coordinate for some particles in a cluster with gaussian energy density 
profile and total mass M « 0.57 as a function of time. Further we can see that the 
particles retain roughly the same radial coordinate for a time longer than 500M (2 com- 
plete orbits). Also, the motion appears to be stable since the particles oscillate radially 
with roughly constant amplitude (see inset). We found that for clusters with gaussian 
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Figure 4.2: Radial coordinates of 7 particles as a function of time for a static cluster with 
gaussian energy density profile. 

energy densities and initial maximum 6M/r = 1.01, the cluster becomes unstable and 
formed a black hole, whereas for clusters with initial maximum 6M/r = 0.97 the cluster 
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was stableQ. This result agrees with the result that Einstein clusters become unstable for 
6M/r = 1 (see || and references therein). 

We have also simulated distorted static clusters for which we specify (1 — e) times 
the angular momentum fixed by Q4.6Q ; this allows us to control the deviation from the 
ideal static case. We have observed results similar to the case just described where the 
perturbation is effectively introduced by the numerical approximation itself; i.e. our 
distorted clusters tend to be stable when 

max (6M(t, r)/r) < 1. (4.7) 

4.3 Numerical Errors. 

There are two principal sources of numerical errors in our code: 1) sampling of the initial 
distribution function with a finite number of particles (statistical error), and 2) approxi- 
mation of the derivatives in the field equations by finite-difference operators (truncation 
error) . 

In the continuum limit — where the number of particles, N, approaches infinity, and 
the mesh spacing, h, goes to — we need to check that these numerical errors go to zero 
to establish that, in this limit, we are solving the continuum system of equations (Vlasov 
equation coupled with Einstein's equations). 

First let us study how the statistical error varies with the number of particles. We 

can estimate the statistical error in one function by subtracting the functions computed 

using two different initial sets of particles (two different sets, each sampled using a Monte 

Carlo method and the same distribution function). Here, the specific function which we 

monitor is 2M(t,r)/r, and we take the L2 norm of the difference to be an estimate of 

the statistical error at each time. In other words, if ^ (2M(t,r)/r) and *- 2 ^ (2M(t,r)/r) 
1 This cluster actually achieves gets 6M/r > 1 during certain parts of the evolution. 
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are calculated from two different sets of particles (generated from the same distribution 
function), then the statistical error estimate at a given time t is computed by: 



As(t , r) = W (2M(t c , r)/r) - ^ (2M(t , r)/r) 
Furthermore, we take the L2 norm of this function: 

As(t ) 



(4.8) 



» 52As(t ,n) 2 /N r 
\ i=l 



(4.9) 



In Figure |4.3| we have plotted these estimates as a function of the number of particles 
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-4.5 



Figure 4.3: Statistical Error. 

at the initial time t = and for fixed mesh spacing, h. We have displayed 9 different 
estimates for the statistical error, and the shaded area of the graph shows bounds for the 
ensemble. The average slope is m = —0.44 ± 0.08 where the uncertainty is the standard 
deviation calculated from the individual slopes. We see that these results roughly agree 
with the expectation that the statistical error should scale as N~ l l 2 . 

Let us turn our attention now to the truncation error. We assume a Richardson 
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Figure 4.4: Truncation (dotted lines) and statistical error estimates (dashed lines) at 
t=0. 



expansion for a mesh function q h (t,r), calculated with a given mesh spacing h: 



q h (t, r) = q(t, r) + Aq s (t, r) + hPe p {t, r) + 0{h p+1 ) 



(4.10) 



where Aq s is the statistical error. We estimate the truncation error in q h (t,r) via 



q h (t,r)-q h '(t,r) 



(4.11) 



where h! < h, specifically we have used log(Zi') ~ —2.4). In this way we use the function 
calculated with hi as a better approximation (a reference) to the analytic solution than 
the one with larger h. It is important that the statistical error Aq s <C h p e p otherwise this 
recipe to calculate the estimate for the truncation error is meaningless since q h (t,r) — 
q h '(t,r) would be dominated by the statistical error. 



In Figure |4.4j we show the base ten logarithm of the norm of the estimates of both 
the truncation and the statistical errors at t — 0. We have plotted the truncation error 
estimates with dotted lines, and the statistical ones with dashed lines. 
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In these calculations we have not fixed the total number of particles, but rather the 
number of particles per grid cell (the value we use in this particular case is 500 particles 
per grid cell on average). The reason for fixing the number of particles per cell, rather 
than the total number of particles, is that the statistical error scales with h (for the 
total number of particles fixed) roughly as h~ l . Increasing the number of particles as 
we decrease h, we compensate for this bias by maintaining roughly constant and small 
statistical error. 

In the part of the Figure where the statistical error is much smaller than the truncation 
error we have calculated the average slope for the truncation error (the solid line in the 
Figure). The resulting slope is m = 2.02 ± 0.08 where again the uncertainty is the 
standard deviation. At t = we thus observe second order convergence which is in 
agrement with the scheme that we have used to compute the geometric quantities. 




Log(h) 

Figure 4.5: Truncation (dotted lines) and statistical error estimates (dashed lines) at 
t=20. 

Figure [4.5| shows the results of our error analysis at t = 20. We see that in this case 
the slope is not 2 but smaller, m = 1.36 ± 0.08. This is consistent with the observation 
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made in section ^T3] that the numerical scheme that we have used for the evolution is only 
first order accurate in h. 



Chapter 5 



Results 



In this chapter we summarize the main results we have obtained. In order to find the 
critical solutions (solutions sitting at the threshold of black hole formation) we have 
performed bisection searches using the total rest mass (sum of the rest masses of the 
individual particles) of the cluster, M Q , as the search parameter p. In section [5.1| we 
discuss some characteristics of the critical solutions. In 5.2 we give some indications 
for the existence of a scaling laws, r ~ — aTn |p — p*\, and we discuss the possibility of 



universality in the model in section [5.3| . The chapter ends with some conclusions and 
some ideas for future work. 



5.1 Is the critical solution Static ? 

In this section we present evidence that, for initial data with non-zero angular momentum, 
as we approach criticality the solution for the space time approaches a static solution. 

We will first focus on one family of initial conditions and afterwards we will try to 
explain what we observed for different families. The initial distributions on equation 
( ggg ) are: 

R{r) oc r 2 e(- (r - ro)2/A ' 2 ) 0(r) (5.1) 

P( Pr ) OC e (-(Pr-Pro) 2 /A pr 2) (5 ^ 2) 

L(l) oc e(- (/ - io)2/A ' 2 ) 0(/) (5.3) 
where p r = p r /a, and 9 is the step or Heaviside function. For r Q = 5, A r = 1, p ro = 0, 
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Pr 



2, l Q = 12, A; = 2 (we will call this family "almost" time symmetric because 



the gaussian for the radial momentum is centered around p r = 0) we found the critical 
solution for a cluster rest mass about M m 1.3. In Figure [5J] we show a few snapshots 
from the evolution of da/dt resulting from initial data which is close to criticality but 
which eventually disperses. At early times, da/dt oscillates, but between t = 117 to 
t = 195 it seems to get close to zero. In the last snapshot (t = 234) we observe how da/dt 
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Figure 5.1: Snapshots of da/dt. 



has become negative, corresponding to dispersal of the particles. We observe similar 
behavior for the time derivatives of all the metric coefficients. 



In order to see how close to zero da/dt becomes, we show in Figure |5.2j , a detail of 
a(t,r) at t = 156 for three different sets of particles. The statistical error is roughly of 
the same order as the function itself (this is not true from r ~ 6 to r ~ 8 where the 
three ensembles agree on a non-zero value for the derivative, however we suspect that 
this feature will decrease if we tune closer to the critical solution and if we use greater 
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resolution). This shows that the derivative of the function could be zero as we anticipate 
but it doesn't give us a definite answer. The value of da/dt could be smaller that the 
statistical error. In order to definitively answer this question we will need to have better 



resolution and decreased statistical errors. As we have commented in section 2.4 an 



X: 




10 

r 



t = 156 

Ensemble 1 

Ensemble 2 
Ensemble 3 



Figure 5.2: Three ensembles of da/dt at t = 156. 



interesting function we monitor is 2M(t,r)/r. This function is a scale invariant quantity 
which tells us a how close we are to trapped surface formation. We plot 

2M(t,r) 



max • 



(5.4) 



in Figure |5.1] . In the inset we show a detail for the period of time when the maximum 
seems to have a roughly constant value. Again, between t = 140 and t = 190 the 
statistical errors are of the same order as the fluctuations in 2M(t,r)/r 

If we accept that the metric coefficients are independent of our time coordinate then 
we have a stationary space-time. If, in addition, the vector N a = d/dt introduced in 
Chapter (Q) is orthogonal to the spatial hypersurfaces, then the space-time is static. 
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Figure 5.3: max r (2M(t,r)/r) as a function of time. 

This then implies that Q a = N a , or that the shift vector, and in particular the shift 
function /3(t,r), is zero. 

In Figure ^| we show the evolution of the shift function. During the period when 
the time derivatives of the metric coefficients are close to zero, the shift function /3(t,r) 
also is close to zero, or at least of the same order as the statistical fluctuations as before. 
Thus the solution for the space-time sitting at the threshold of black hole formation 
could be a static solution (if it is not static, the amplitude of the derivatives and the shift 
function have small amplitudes relative to the statistical error). We also note that the 
total current density j r also tends to be of the order of the statistical fluctuations, but 
that the component of the stress energy tensor is not zero in the critical regime, as 
we can see in Figure |5.5| . Therefore, on average there are the same number of particles 
with positive (outwards) and negative (inwards) r components of the 4-momentum but 
the absolute value of p r is not zero on average. 
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Figure 5.4: Shift function as a function of time. 

This result and the fact that the maximum value of 2M(r)/r for this critical solution 
is ~ 0.76 indicate that this cannotbe one of the Einstein clusters we have considered pre- 
viously since there are no equilibrium clusters (either stable or unstable) with maximum 
2M(r)/r larger than 2/3. 
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Figure 5.5: Evolution of S r (t,r). 



5.2 Time Scaling 



A Type I critical solution usually approaches a static or a periodic solution. As was 
argued in the previous section, we believe the critical solution in the current case is 
static. As we approach p* the dynamical solution spends more and more time close to 
this putative static solution. 

We have calculated the time t that each of the dynamical solutions spends inside some 
radius r = r Q (in particular r = 6). Figure |5.6| shows the plot of r = t — t c versus ln(p) 
(natural logarithm) where t is the time that the solution with search parameter p spends 
inside r = r Q , and t c is the same quantity for the solution closest to criticallity. 

We show three different ensembles of particles for the gaussian family given by equa- 
tions ( |5.1[ - |5.3| ) and again using "almost time symmetric" initial data. We also have nor- 
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Figure 5.6: Scaling law for the time. The error for each value of a is the standard 
deviation of the slope computed using a least squares fit. 

malized the critical solution in such a way that it has ADM mass equal to 1 at infinity^, 
i.e.: 

r -> r/M c (f\r max ) (5.5) 
t _> t/M c (r,r max ), (5.6) 

where M c (t*,r max ) is the value of the mass aspect function at r = r max for the solution 
closest to criticality during the critical regime. We can see that there is a rough linear 
relation between the time r and ln|p— p*\ with slope —(5.2 ±0.2) (the error is an estimate 
of the statistical error between different sets of particles): 

r (5.2 ±0.2) In \p - p*\ (5.7) 

Qualitatively this agrees with other type I critical solutions. 
lr rhis will be useful when we compare with other families of initial data. 
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5.3 Different Families of Initial Data 

The results we showed in the previous section are for Gaussian initial data given by 



equations ||5.1[ - 15.3|| with p ro = ( "almost time symmetric" initial data) and l Q = 7. In 
this section we show that similar results are found for other families of initial data. We 
have chosen initial data with p ro = —4 and l Q — 3, l — 5, l Q = 7, l Q = 12, (with same 
r D = 5, A r = 1, Ap r = 2, Ai = 2 as before) and find that the critical solution for each 
of these cases appears to approach a static solution. For smaller values of the angular 
momentum the mass in the critical solution gets concentrated closer to the origin, which 
makes it more difficult to resolve the solution with a uniform finite-difference grid. We 
also have evolved an initial family with the following one particle distributions: 

3(r) (5.8) 
V)) (5-9) 



R(r) 


OC (l 


— tanh ((r — 


P(Pr) 


oc (1 


— tanh ((p r - 


L(l) 


oc (l 


— tanh ((/ — 


TO — 


-4 A- 


= 2, l = 7, 



/A,) 2 ) 9(0 (5.10) 

with r Q = 5,A r = 1, p ro = —4, A pr — 2, l — 7, A^ = 2. For this distribution we also 
find that the solution with small p — p* (p being the total rest mass of the cluster as 
before) spends some time close to an apparently static solution. In Figure 1577] we show 



profiles for 2M(r)/r for the different families, each separate profile being selected from 
the corresponding period of near-critical evolution. Since different initial conditions set 
different scales for the problem we have normalized the results such that the total ADM 
mass of the solutions to which they approach is 1 at infinity (rescaling given by equations 
( |5.5| - |5.6|) ). We can see that after normalization, the peak of 2M(r)/r is roughly at the 



same r* = 2.3. We can also see that the better a solution is resolved, the closer it conforms 
to the best resolved solution (the Gaussian with l = 7 and "almost time symmetric" 
initial data). This is an indication that the critical solution might be universal although, 
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we again need better resolution and more particles to be able to give a definitive answer. 

We are also interested in seeing the importance of the distribution of angular momen- 
tum on the critical solutions. In Figure |5.8| we show r 2 S a a (r) where S a a (r) is the function 
defined by equations ( |2.58|) and (|2.56|) for the different families of initial data during the 
critical regime. This function, r 2 S a a (r), is a dimensionless quantity which measures the 
square of the angular momentum of the distribution of particles. We have rescaled the 
radial coordinate (and time) in the same way as in Figure |5.7|. We can see that there is 
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Figure 5.8: r S a a (r) for the different families during the critical regime. 



no apparent agreement between the functions calculated from different families of initial 
data. More work needs to be done in order to see what is the dependence of the critical 
solution with respect the angular momentum distribution. Moreover, estimates for the 
truncation and statistical errors for each calculation of the critical solutions must be 



Chapter 5. Results 



49 



calculated to properly compare the solutions. 
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Figure 5.9: Scaling behaviour for different families of initial data . 

We have also estimated a defined by 

r = — a\n\p — p*\ 



(5.11) 



for the different families. In Figure |5.9| we show the scaling behaviour for the different 
families of initial data that we have studied in this thesis. The error for each value of a 
in the Figure is, as beforre, the standard deviation of the slope computed using a least 



squares fit. We have collected the values that we have obtained for a in Table [57L 



5.4 Conclusions and Future Work 



We have studied critical behavior at the threshold of black hole formation for collisionless 
matter. Our results indicate that, for families with non-zero angular momentua, the 
critical solution could be static with non-zero radial particle momenta. We have found 



Chapter 5. Results 



50 



Family 


G 


Gaussian, l Q — 7 (set 1) 


5.0 ± 0.2 


Gaussian, l = 7 (set 2) 


5.0 ± 0.2 


Gaussian, Z D = 3 


5.7 ±0.2 


Gaussian, l — 5 


5.5 ±0.2 


Gaussian, Z Q = 12 


4.9 ±0.2 


Gaussian, Z Q = 12 t.s. (set 1) 


5.1 ±0.2 


Gaussian, l Q = 12 t.s. (set 2) 


5.3 ±0.2 


Gaussian, Z Q = 12 t.s. (set 3) 


5.2 ±0.2 


Tanh, l Q = 7 


5.9 ±0.2 



Table 5.1: Values of time scaling exponent for the different families (t.s. stands for 
"almost time symmetric"). The errors above are assumed to be the same as the statistical 
error for the t.s. family (see equation (|5.7|) ). 

evidence for a life-time scaling law which is to be expected for Type I critical solutions. In 
order to answer all these questions more rigorously we would need some way to increase 
our resolution where it is needed. Some kind of adaptive code (such as the ones used in 
finite difference studies [|1|) would be useful. 

It may be that the development of a finite difference code to solve the Vlasov equation 
directly in phase space would be the best way to attack this problem. This would allow us 
to incorporate the adaptivity that we need, as well as providing us with better-understood 
convergence properties. 
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Appendix A 



Summary of Equations in Maximal- Areal Coordinates 



A.l Field equations 



where 



and 



a! I -a 2 3 



+ -ra 2 Kf + Airarp 



a 2r 2 

K f = J_ K e 

r a 

a" = a '(---) + ^(a 2 -l + 2r-) + Anaa^ + S e e - 3p) 



a r r 2 a a 2 



(3 = arK 6 e 



1 N 0* 

i N i 2 

i N v . 



W(x - Xi) 



4nAr ~^ r, 

l — \x — Xi\/h : \x — Xi\ < h 
: otherwise 
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A. 2 Evolution Equations 



dp r _ da t dp a dap 2 l 2 a 
dt or or a 6 or p l p t r i 



dt r (a 2 p* ^) ^ ^ tanOp 1 (r 4 ^ ) ^ ^ 

dO ap e 



dt p 



(A.25) 



P* 2 = - 1 T-A l --P e2 ) (A-26) 
sm 6 \r 4 J 

dd> ap^ , . 

-r = -zr A.27 

dt f y 1 



Appendix B 

Direct Approach To Solutions of Einstein- Vlasov System 



B.l Field Equations 

The distribution function is defined by: 

f(t,r,f,l) = ^ F , (B.28) 

where N is the number of particles and V is the volume in phase space. Then the stress 
energy tensor can be calculated via: 

T^= [ p»p»f(t,r,p r ,l)-dV p (B.29) 
J m 

where the volume element in momentum space, V p , is restricted by 

- m 2 = pV (B.30) 
to the following expression (see [fD|): 

dp r dp dp^ 

dV„ = m — - — (B.31) 

We can introduce momentum coordinates adapted to the symmetry: 

I 2 = P 2 e + A- ft (B.32) 
sin V 

1 pe sin 

ip = tan- 1 - B.33 

and p* given, as before, by: 



pt = J m * + f 2 +^ (B.34) 
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The volume element dV p in these variables is: 



dp r dl 2 dip 

= m ^^T (B.35) 



The stress-energy tensor components are 

f —t £ j „2 , fa 



T t = ^ f P'fdPrdl 2 + I p r fdp r dl 2 (B.36) 

/ prfdp r dl 2 (B.37) 

2 

/ p r /^ r rf/ 2 + / ^fdp r dl 2 (B.38) 



rpt _ 71 



a?r 2 ot 



rpr = ~fa 



aa 2 r 2 J ~ a^r* J p 



T e= T t = 2^/^ 2 (R39) 



We are able then to construct the source terms: 



P = —2 I f(t,r,p r ,l)Jm 2 + ^ + ^dp r dl 2 (B.40) 
ar z J V a r 2 

j r = — r / f{t,r,p r ,l)p r dp r dl 2 (B.41) 
ar z J 

SI = JL/ - tif(t,r,Pr,l) 2 (R42) 

S e e = ^ / ; ££g^J rf/ 2 ^ (B.43) 



ar 4 



^m 2 + p 2 /a 2 + l 2 /r 



B.2 Vlasov Equation 



If the system is collisionless, we have Liouville's theorem concerning the conservation 
of the volume in phase space, V, and we also know that the particle number N is a 
conserved quantity. Therefore we can write the conservation of the distribution function, 

/ = N/V, as: 

f - = (B.44) 
dr 
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This equation is known as the "collisionless Boltzmann's equation" or the "kinetic equa- 
tion". Expanding the total derivative we have: 

dx a df dp a df 



+ 







^B.45) 



dr dx a Or dp c 

where x a are the space-time variables, and p a are the momentum variables. In spherical 
symmetry this equation reduces to: 

df dr df dp r df 



dt ~^ dt dr ~^ dt dp r 



(B.46) 



where dr/dt and dp r /dt are given by equations ( 2.77|) and Q2.78 ). This gives: 





df_ , apr 
dt 



df da_ t dp adap r 2 al 2 df 
a^p 1 - dr dr dr er dr p 1 p^r 6 dp r 



(B.47) 



Appendix C 



Summary of Equations in Polar-Areal Gauge Coordinates 



In this section we present the equations for a second choice of coordinates. Specifically 
we again adopt areal radial coordinates (implying that the proper surface area of the 
spheres centered at r = and with radius r will be 4nr 2 ), but fix the time coordinate 
using polar slicing. This slicing condition is implemented by demanding: 

TrK = K\ = K r r (C.48) 

which implies K% = = K\ = = 0. Using expression fl2.39|) we can express 



K\ as: 



K\ = —((3(rb)'-rb) = 0, (C.49) 
rba 



we have: 

(3{rb)' = rb (C.50) 
Choosing areal coordinates, b = 1, 6 = 0, the above equation implies: 

(3 = 0. (C.51) 

The resulting metric is: 

ds 2 = -a 2 (t, r)dt 2 + a 2 (t, r)dr 2 + r 2 dfl (C.52) 
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we then have the following geometrical quantities: 



R = 2 



-2(l/a)' + a/r(l- l/a 2 )l /(ar) 



i? r r = 2a'/ (ra 3 ) 



J R% = i?^ = l/(ar 2 ) (a-(r/a)') 



K r r = K = —a/ (aa) 



K% = K^ = 



(C.53) 



K % r\i = — (a/ (aa))' — 2a/ (race) K\ r = K r rr = — (a/ (aa)) 

a\ r \ r = 1/a (a'/a)' a^ e = = a'/ (a 2 r) 

Briefly we summarize the resulting equations of motion: 
Momentum Constraint: 

a = —Anraajr (C.54) 

Hamiltonian Constraint: 

a 1 = 4irra 3 p-— (a 2 - 1) (C.55) 
2r 

Evolution of the 3-metric 

= -2ag ik K* (C.56) 

Evolution of the extrinsic curvature 

K\ = -- (-) +^ + 4ira(S-p)-8naS r r (C.57) 
a \ a J rar 

K% = = -^ + ^(a--+S) + 4ira(S-p)-8naS e (C.58) 
a 2 r ar" 1 \ a a 2 J 

K% = = -^- + ^-( a --+r^-) +4ira(S-p)-87raS% (C.59) 
a 2 r ar 1 \ a a 1 J 

where the two last equations are zero by the choice of slicing. 



Appendix D 



Pseudo-code 



routine vlasov 

# ri := position of the i-th particle 

# \Pfj]i := momentum of i-th the particle 

# [T^jj := T/f component on the j grid point 

# aj, «j and (3j are the metric coefficients 

# Mj := mass aspect function 

# Generate the initial positions r,- L and velocities [p M ]j for the particles 
read parameters of the distribution function 

for % 1 ... N particles 

Generate T{ and \pp]i using Monte Carlo for a given distribution function 
end for 
t = 

# Time loop 
while t < t max 

for j ^-■■■-^gridpoints 

Calculate [Ttf]j from and \p^\i 
end for 

Solve for ctj, (3j and ctj on the mesh 
# Check apparent horizon formation 
if 2Mj /rj = 1 stop 
Ouput ctj, Pj and aj 
for i 1 ... Np ar n c i es 

Calculate [r^]j at i = i + 1/2 (it from a,-, (3j and ay 
end for 
t = t + dt 
for i 1 ... Np ar n c i es 

Integrate geodesic equations to find new and \pp\i 

if ri > r max Kill i-th particle, N partides = N partides - 1 

if rj < Reflect the particle 

if Npartides = Stop 

end for 
end while 
end routine 
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